library(tidyverse)
library(openxlsx)
library(ggplot2)
library(ggpubr)
library(gridExtra)
library(survival)
library(simPH)
library(texreg)
library(modelsummary)
library(kableExtra)
library(tinytex)
library(ggallin)
library(RColorBrewer)
library(ggrepel)
library(dplyr)

##################################################################################
####################### Data files ##############################################
# Share of independent MPs
ind <- read.xlsx("IndMPShareCombined.xlsx", detectDates = TRUE)

# Main dataset
new <- readRDS("ReplicationData.rds") 

# Codebook
# country - country (RO1: RO PR period; RO2: RO MMP period)
# term - parliamentary term
# country_term - country+parl. term
# name - MP name
# date - date
# start and stop - number of days since becoming independent
# spellid - ID of the spell of independence
# affilevent - did MP end their spell of independence by joining an existing or new group?
# event - numerical variable for the end of spell of independence period: 0: no entry or PPG creation; 1: PPG entry; 2: PPG creation
# entryevent - did MP end their spell of independence by joining an existing group?
# PPGcrevent - did MP end their spell of independence by joining a new group?
# lasday - number of days from becoming an independent until entering a group or the end of the tenure in parliament
# prominence - does MP have a strong personal vote?
# MinorParty - was an MP a member of minor party that does not have its own parliamentary group?
# govmaj - the share of MPs in government parties
# logterms - log of the number of parliamentary terms MP has served
# dist2 - policy proximity between the MP and two closest parl. parties
# collapse - did the former PPG of the MP collapse?
# priorel - number of days since the prior election
# nextel - number of days until the next election
# Expulsion - was MP expelled from their former group?
# coswitcher - share of co-leaver MPs from the total number of MPs
# relcoswitcher - share of co-leaver MPs from the total number of MPs minus the share of MPs required to establish a PPG
# eldate - date of the prior election
# eldate_next - date of the next election
# logprefshare - log of the share of preference votes (from the total number of votes cast in a district) by the MP, in the prior election
# govformation - government formation period
# smdsurplus - (share of MPs' vote share in SMD contest / share of party's support in the PR district) - 1
# preel - pre-election period
# lt, pl, ro1, ro2 - country identifiers
# variables with "t" - products of logged number of days as independent and respective substantive variables

lt <- new[new$country=="LT",]
pl <- new[new$country=="PL",]
ro1 <- new[new$country=="RO1",]
ro2 <- new[new$country=="RO2" ,]

# List of independence spells
spells <- new %>% filter(duplicated(spellid)==FALSE) %>% dplyr::select(-event)
tmp <- new %>% dplyr::select(spellid,lastday,event,stop) %>% filter(lastday==stop) %>% dplyr::select(-stop)
spells <- left_join(spells,tmp)
##################################################################################
##################################################################################


##################################################################################
################################ FIGURE 1 #######################################
ind$ind_share <- ind$ind_share*100
ggplot(ind,aes(x=day,y=ind_share)) + 
  geom_line(na.rm = TRUE) + 
  theme_bw(base_family = "serif") +
  theme(plot.caption = element_text(hjust = 0)) +
  facet_grid(country~.) + 
  labs(x="", y="Share of independent MPs", caption = "Source: Own elaboration")
ggsave("figure_1.png", height = 5, width = 7, dpi = 300)
##################################################################################
##################################################################################


##################################################################################
################################# FIGURE 2 ######################################
temp <- new %>% filter(lastday==stop & affilevent==1) %>%
  mutate(event=ifelse(PPGcrevent==1,"PPG creation","Entry"))
templt <- temp %>% filter(country=="LT")
temppl <- temp %>% filter(country=="PL")
tempro1 <- temp %>% filter(country=="RO1")
tempro2 <- temp %>% filter(country=="RO2")

noswitch <- new %>% filter(lastday==stop & affilevent==0)
noswitchlt <- nrow(noswitch[noswitch$country=="LT",])
noswitchpl <- nrow(noswitch[noswitch$country=="PL",])
noswitchro1 <- nrow(noswitch[noswitch$country=="RO1",])
noswitchro2 <- nrow(noswitch[noswitch$country=="RO2",])

data <- data.frame(Country = c("Lithuania", "Poland", "Romania-PR", "Romania-MMP"),
                   PPG_entry = c(table(templt$event)[1], table(temppl$event)[1],   
                                 table(tempro1$event)[1],table(tempro2$event)[1]),
                   PPG_creation = c(table(templt$event)[2], table(temppl$event)[2],   
                                    table(tempro2$event)[2], table(tempro2$event)[2]),
                   No_Event = c(noswitchlt,noswitchpl,noswitchro1,noswitchro2)) %>%
  rename_with( ~str_replace(., "_", " ")) %>%
  pivot_longer(cols=c(2:4), names_to="ppg_type", values_to="ppg_freq")

data <- data %>% 
  group_by(Country) %>%
  mutate(
    total = sum(ppg_freq),
    percentage = ppg_freq / total * 100,
    label = paste0(ppg_freq, " (", round(percentage), "%)")
  ) %>%
  ungroup()

# Create the plot showing proportions with count and percentage labels
freq <- ggplot(data, aes(x = Country, y = ppg_freq, fill = ppg_type)) +
  geom_col(position = position_fill()) +
  scale_fill_manual(values = c("darkgray", "steelblue", "tomato3")) +
  theme_bw(base_family = "serif") +
  theme(legend.title = element_blank()) +
  labs( y="Share of events", x="Country") +
  geom_text(aes(label = label), 
            position = position_fill(vjust = 0.5), 
            color = "white", fontface = "bold", size = 3) +
  # Format y-axis as percentages
  scale_y_continuous(labels = scales::percent_format())

ggsave("figure_2.png", height = 4, width = 6, dpi = 300)
##################################################################################
##################################################################################


##################################################################################
################################# FIGURE 3 #######################################
# Figure 3: Time until entering or creating a group
temp <- rbind(cbind(templt[c("stop","event")],country="LT"),
              cbind(temppl[c("stop","event")],country="PL"),
              cbind(tempro1[c("stop","event")],country="RO-PR"),
              cbind(tempro2[c("stop","event")],country="RO-MMP"))
ggplot(temp, aes(x = stop)) +
  geom_histogram(binwidth=5,fill="white",color="black") +
  facet_grid(event~country, scales = "free") +
  theme_bw(base_family = "serif") +
  theme(plot.caption = element_text(hjust = 0)) +
  labs(title="", y="Count", x="Days as independent MP", caption = "Source: Own elaboration")
ggsave("figure_3.png", height = 5, width = 8, dpi = 700)
##################################################################################
##################################################################################



##################################################################################
###################### Models and Figures 4 and 6 ################################

##################################################################################
#### Regression models ###########
# without distance variable, ENTRY
m1lt <- coxph(Surv(start,stop,entryevent) ~ prominence + MinorParty + govmaj + collapse + Expulsion + coswitcher + preel+ govformation,lt,id=spellid,cluster=spellid, robust=T)
cox.zph(m1lt,transform="identity", global=F)
m1pl <- update(m1lt,data=pl)
cox.zph(m1pl,transform="identity", global=F)
m1pl <- update(m1pl,.~.+tMinorParty+tgovmaj)
simMinorParty <- coxsimtvc(m1pl, "MinorParty", "tMinorParty",qi = "Relative Hazard", from=1, to=365, Xj=1,tfun="log",nsim=1000)
simGG(simMinorParty,xlab="Days",legend=FALSE)
ggsave("figure_3_Appendix.png", height = 6, width = 7, dpi = 300)
m1ro1 <- update(m1lt,data=ro1)
cox.zph(m1ro1,transform="identity", global=F)
m1ro1 <- update(m1ro1,.~.+tcollapse+tcoswitcher+tpreel+tgovformation)
m1ro2 <- update(m1lt,data=ro2)
cox.zph(m1ro2,transform="identity", global=F)

# without distance variable, PPG CREATION
m3lt <- coxph(Surv(start,stop,PPGcrevent) ~ prominence + MinorParty + govmaj + collapse + Expulsion + relcoswitcher + preel + govformation,lt,id=spellid,cluster=spellid, robust=T)
cox.zph(m3lt,transform="identity", global=F)
m3pl <- update(m3lt,data=pl)
cox.zph(m3pl,transform="identity", global=F)
m3ro1 <- update(m3lt,data=ro1)
#cox.zph(m3ro1,transform="identity", global=F)
m3ro2 <- update(m3lt,data=ro2)
cox.zph(m3ro2,transform="identity", global=F)

# with distance variable, ENTRY
m2lt <- update(m1lt,.~.+dist2)
cox.zph(m1lt,transform="identity", global=F)
m2pl <- update(m2lt,data=pl)
cox.zph(m2pl,transform="identity", global=F)
m2pl <- update(m2pl,.~.+tMinorParty+tgovmaj+tpreel)
m2ro1 <- update(m2lt,data=ro1)
cox.zph(m2ro1,transform="identity", global=F)
m2ro1 <- update(m2ro1,.~.+tgovmaj+tcollapse+trelcoswitcher)
m2ro2 <- update(m2lt,data=ro2)
cox.zph(m2ro2,transform="identity", global=F)

# with distance variable, ENTRY
m4lt <- update(m3lt,.~.+dist2)
cox.zph(m4lt,transform="identity", global=F)
m4pl <- update(m4lt,data=pl)
cox.zph(m4pl,transform="identity", global=F)
m4ro1 <- update(m4lt,data=ro1)
#cox.zph(m4ro1,transform="identity", global=F)
m4ro2 <- update(m4lt,data=ro2)
#cox.zph(m4ro2,transform="identity", global=F)
##################################################################################


##################################################################################
######################## Simulations for Figures 4 and 6 ########################
# (for time-dependent variables, get the value on day 148: median observed value of non-affiliation period)

##################################################################################
# Values of predictor variables for Figures 4 and 6
govmaj_lower <- mean(new$govmaj) - sd(new$govmaj)
govmaj_upper <- mean(new$govmaj) + sd(new$govmaj)
relcoswitcher_lower <- mean(spells$relcoswitcher) - sd(spells$relcoswitcher)
relcoswitcher_upper <- mean(spells$relcoswitcher) + sd(spells$relcoswitcher)
coswitcher_lower <- max(0,mean(spells$coswitcher) - sd(spells$coswitcher))
coswitcher_upper <- mean(spells$coswitcher) + sd(spells$coswitcher)
##################################################################################

##################################################################################
#####################PPG entry as the outcome variable###########################
# Lithuania 
sim_prominence_entry_lt <- cbind(simGG(coxsimLinear(m1lt, b="prominence", Xj=1, Xl=0, nsim=1000, ci=0.95))$data[c("Median","Min","Max")])
sim_MinorParty_entry_lt <- cbind(simGG(coxsimLinear(m1lt, b="MinorParty", Xj=1, Xl=0, nsim=1000, ci=0.95))$data[c("Median","Min","Max")])
sim_govmaj_entry_lt <- cbind(simGG(coxsimLinear(m1lt, b="govmaj", Xj=govmaj_upper, Xl=govmaj_lower, nsim=1000, ci=0.95))$data[c("Median","Min","Max")])
sim_collapse_entry_lt <- cbind(simGG(coxsimLinear(m1lt, b="collapse", Xj=1, Xl=0, nsim=1000, ci=0.95))$data[c("Median","Min","Max")])
sim_Expulsion_entry_lt <- cbind(simGG(coxsimLinear(m1lt, b="Expulsion", Xj=1, Xl=0, nsim=1000, ci=0.95))$data[c("Median","Min","Max")])
sim_relcoswitcher_entry_lt <- cbind(simGG(coxsimLinear(m1lt, b="coswitcher", Xj=coswitcher_upper, Xl=coswitcher_lower, nsim=1000, ci=0.95))$data[c("Median","Min","Max")])
sim_preel_entry_lt <- cbind(simGG(coxsimLinear(m1lt, b="preel", Xj=1, Xl=0, nsim=1000, ci=0.95))$data[c("Median","Min","Max")])
sim_govformation_entry_lt <- cbind(simGG(coxsimLinear(m1lt, b="govformation", Xj=1, Xl=0, nsim=1000, ci=0.95))$data[c("Median","Min","Max")])
sim_dist2_entry_lt <- cbind(simGG(coxsimLinear(m2lt, b="dist2", Xj=1, Xl=0, nsim=1000, ci=0.95))$data[c("Median","Min","Max")])


# Poland
sim_prominence_entry_pl <- cbind(simGG(coxsimLinear(m1pl, b="prominence", Xj=1, Xl=0, nsim=1000, ci=0.95))$data[c("Median","Min","Max")])
sim_MinorParty_entry_pl_full <- simGG(coxsimtvc(m1pl, b="MinorParty", btvc="tMinorParty", Xj=1, 
                                                Xl=0, from=1, to=365, nsim=1000, ci=0.95, tfun="log"))$data
sim_MinorParty_entry_pl_full <- sim_MinorParty_entry_pl_full[order(sim_MinorParty_entry_pl_full$Time),]
sim_MinorParty_entry_pl <- cbind(sim_MinorParty_entry_pl_full[c("Median","Min","Max")][148,])
sim_govmaj_entry_pl_full <- simGG(coxsimtvc(m1pl, b="govmaj", btvc="tgovmaj", Xj=govmaj_upper, 
                                            Xl=0, from=1, to=365, nsim=1000, ci=0.95, tfun="log"))$data
sim_govmaj_entry_pl_full <- sim_govmaj_entry_pl_full[order(sim_govmaj_entry_pl_full$Time),]
sim_govmaj_entry_pl <- cbind(sim_govmaj_entry_pl_full[c("Median","Min","Max")][148,])
sim_collapse_entry_pl <- cbind(simGG(coxsimLinear(m1pl, b="collapse", Xj=1, Xl=0, nsim=1000, ci=0.95))$data[c("Median","Min","Max")])
sim_Expulsion_entry_pl <- cbind(simGG(coxsimLinear(m1pl, b="Expulsion", Xj=1, Xl=0, nsim=1000, ci=0.95))$data[c("Median","Min","Max")])
sim_relcoswitcher_entry_pl <- cbind(simGG(coxsimLinear(m1pl, b="coswitcher", Xj=coswitcher_upper, Xl=coswitcher_lower, nsim=1000, ci=0.95))$data[c("Median","Min","Max")])
sim_preel_entry_pl <- cbind(simGG(coxsimLinear(m1pl, b="preel", Xj=1, Xl=0, nsim=1000, ci=0.95))$data[c("Median","Min","Max")])
sim_govformation_entry_pl <- cbind(simGG(coxsimLinear(m1pl, b="govformation", Xj=1, Xl=0, nsim=1000, ci=0.95))$data[c("Median","Min","Max")])
sim_dist2_entry_pl <- cbind(simGG(coxsimLinear(m2pl, b="dist2", Xj=1, Xl=0, nsim=1000, ci=0.95))$data[c("Median","Min","Max")])


# Romania1
sim_prominence_entry_ro1 <- cbind(simGG(coxsimLinear(m1ro1, b="prominence", Xj=1, Xl=0, nsim=1000, ci=0.95))$data[c("Median","Min","Max")])
sim_MinorParty_entry_ro1 <- cbind(simGG(coxsimLinear(m1ro1, b="MinorParty", Xj=1, Xl=0, nsim=1000, ci=0.95))$data[c("Median","Min","Max")])
sim_govmaj_entry_ro1 <- cbind(simGG(coxsimLinear(m1ro1, b="govmaj", Xj=govmaj_upper, Xl=govmaj_lower, nsim=1000, ci=0.95))$data[c("Median","Min","Max")])
sim_collapse_entry_ro1_full <- simGG(coxsimtvc(m1ro1, b="collapse", btvc="tcollapse", Xj=1,
                                               Xl=0, from=1, to=365, nsim=100, ci=0.95, tfun="log"))$data
sim_collapse_entry_ro1_full <- sim_collapse_entry_ro1_full[order(sim_collapse_entry_ro1_full$Time),]
sim_collapse_entry_ro1 <- cbind(sim_collapse_entry_ro1_full[c("Median","Min","Max")][148,])
sim_Expulsion_entry_ro1 <- cbind(simGG(coxsimLinear(m1ro1, b="Expulsion", Xj=1, Xl=0, nsim=1000, ci=0.95))$data[c("Median","Min","Max")])
sim_relcoswitcher_entry_ro1_full <- simGG(coxsimtvc(m1ro1, b="coswitcher", btvc="tcoswitcher", 
                                                    Xj=coswitcher_upper, Xl=coswitcher_lower, from=1, to=365, 
                                                    nsim=1000, ci=0.95, tfun="log"))$data
sim_relcoswitcher_entry_ro1_full <- sim_relcoswitcher_entry_ro1_full[order(sim_relcoswitcher_entry_ro1_full$Time),]
sim_relcoswitcher_entry_ro1 <- cbind(sim_relcoswitcher_entry_ro1_full[c("Median","Min","Max")][148,])
sim_preel_entry_ro1_full <- simGG(coxsimtvc(m1ro1, b="preel", btvc="tpreel", Xj=1,
                                            Xl=0, from=1, to=365, nsim=1000, ci=0.95, tfun="log"))$data
sim_preel_entry_ro1_full <- sim_preel_entry_ro1_full[order(sim_preel_entry_ro1_full$Time),]
sim_preel_entry_ro1 <- cbind(sim_preel_entry_ro1_full[c("Median","Min","Max")][148,])
sim_govformation_entry_ro1_full <- simGG(coxsimtvc(m1ro1, b="govformation", btvc="tgovformation", Xj=1,
                                            Xl=0, from=1, to=365, nsim=1000, ci=0.95, tfun="log"))$data
sim_govformation_entry_ro1_full <- sim_govformation_entry_ro1_full[order(sim_govformation_entry_ro1_full$Time),]
sim_govformation_entry_ro1<- cbind(sim_govformation_entry_ro1_full[c("Median","Min","Max")][148,])
sim_govformation_entry_ro1 <- cbind(simGG(coxsimLinear(m1ro1, b="govformation", Xj=1, Xl=0, nsim=1000, ci=0.95))$data[c("Median","Min","Max")])
sim_dist2_entry_ro1 <- cbind(simGG(coxsimLinear(m2ro1, b="dist2", Xj=1, Xl=0, nsim=1000, ci=0.95))$data[c("Median","Min","Max")])


# Romania2
sim_prominence_entry_ro2 <- cbind(simGG(coxsimLinear(m1ro2, b="prominence", Xj=1, Xl=0, nsim=1000, ci=0.95))$data[c("Median","Min","Max")])
sim_MinorParty_entry_ro2 <- cbind(simGG(coxsimLinear(m1ro2, b="MinorParty", Xj=1, Xl=0, nsim=1000, ci=0.95))$data[c("Median","Min","Max")])
sim_govmaj_entry_ro2 <- cbind(simGG(coxsimLinear(m1ro2, b="govmaj", Xj=govmaj_upper, Xl=govmaj_lower, nsim=1000, ci=0.95))$data[c("Median","Min","Max")])
sim_collapse_entry_ro2 <- cbind(simGG(coxsimLinear(m1ro2, b="collapse", Xj=1, Xl=0, nsim=1000, ci=0.95))$data[c("Median","Min","Max")])
sim_Expulsion_entry_ro2 <- cbind(simGG(coxsimLinear(m1ro2, b="Expulsion", Xj=1, Xl=0, nsim=1000, ci=0.95))$data[c("Median","Min","Max")])
sim_relcoswitcher_entry_ro2 <- cbind(simGG(coxsimLinear(m1ro2, b="coswitcher", Xj=coswitcher_upper, Xl=coswitcher_lower, nsim=1000, ci=0.95))$data[c("Median","Min","Max")])
sim_preel_entry_ro2 <- cbind(simGG(coxsimLinear(m1ro2, b="preel", Xj=1, Xl=0, nsim=1000, ci=0.95))$data[c("Median","Min","Max")])
sim_govformation_entry_ro2 <- cbind(simGG(coxsimLinear(m1ro2, b="govformation", Xj=1, Xl=0, nsim=1000, ci=0.95))$data[c("Median","Min","Max")])
sim_dist2_entry_ro2 <- cbind(simGG(coxsimLinear(m2ro2, b="dist2", Xj=1, Xl=0, nsim=1000, ci=0.95))$data[c("Median","Min","Max")])
##################################################################################



##################################################################################
#####################PPG creation as the outcome variable#########################

# Lithuania
sim_prominence_ppg_lt <- cbind(simGG(coxsimLinear(m3lt, b="prominence", Xj=1, Xl=0, nsim=1000, ci=0.95))$data[c("Median","Min","Max")])
sim_MinorParty_ppg_lt <- cbind(simGG(coxsimLinear(m3lt, b="MinorParty", Xj=1, Xl=0, nsim=1000, ci=0.95))$data[c("Median","Min","Max")])
sim_govmaj_ppg_lt <- cbind(simGG(coxsimLinear(m3lt, b="govmaj", Xj=govmaj_upper, Xl=govmaj_lower, nsim=1000, ci=0.95))$data[c("Median","Min","Max")])
sim_collapse_ppg_lt <- cbind(simGG(coxsimLinear(m3lt, b="collapse", Xj=1, Xl=0, nsim=1000, ci=0.95))$data[c("Median","Min","Max")])
sim_Expulsion_ppg_lt <- cbind(simGG(coxsimLinear(m3lt, b="Expulsion", Xj=1, Xl=0, nsim=1000, ci=0.95))$data[c("Median","Min","Max")])
sim_relcoswitcher_ppg_lt <- cbind(simGG(coxsimLinear(m3lt, b="relcoswitcher", Xj=relcoswitcher_upper, Xl=0, nsim=1000, ci=0.95))$data[c("Median","Min","Max")])
sim_preel_ppg_lt <- cbind(simGG(coxsimLinear(m3lt, b="preel", Xj=1, Xl=0, nsim=1000, ci=0.95))$data[c("Median","Min","Max")])
sim_govformation_ppg_lt <- cbind(simGG(coxsimLinear(m3lt, b="govformation", Xj=1, Xl=0, nsim=1000, ci=0.95))$data[c("Median","Min","Max")])
sim_dist2_ppg_lt <- cbind(simGG(coxsimLinear(m4lt, b="dist2", Xj=1, Xl=0, nsim=1000, ci=0.95))$data[c("Median","Min","Max")])


# Poland
sim_prominence_ppg_pl <- cbind(simGG(coxsimLinear(m3pl, b="prominence", Xj=1, Xl=0, nsim=1000, ci=0.95))$data[c("Median","Min","Max")])
sim_MinorParty_ppg_pl <- cbind(simGG(coxsimLinear(m3pl, b="MinorParty", Xj=1, Xl=0, nsim=1000, ci=0.95))$data[c("Median","Min","Max")])
sim_govmaj_ppg_pl <- cbind(simGG(coxsimLinear(m3pl, b="govmaj", Xj=govmaj_upper, Xl=govmaj_lower, nsim=1000, ci=0.95))$data[c("Median","Min","Max")])
sim_collapse_ppg_pl <- cbind(simGG(coxsimLinear(m3pl, b="collapse", Xj=1, Xl=0, nsim=1000, ci=0.95))$data[c("Median","Min","Max")])
sim_Expulsion_ppg_pl <- cbind(simGG(coxsimLinear(m3pl, b="Expulsion", Xj=1, Xl=0, nsim=1000, ci=0.95))$data[c("Median","Min","Max")])
sim_relcoswitcher_ppg_pl <- cbind(simGG(coxsimLinear(m3pl, b="relcoswitcher", Xj=relcoswitcher_upper, Xl=0, nsim=1000, ci=0.95))$data[c("Median","Min","Max")])
sim_preel_ppg_pl <- cbind(simGG(coxsimLinear(m3pl, b="preel", Xj=1, Xl=0, nsim=1000, ci=0.95))$data[c("Median","Min","Max")])
sim_govformation_ppg_pl <- cbind(simGG(coxsimLinear(m3pl, b="govformation", Xj=1, Xl=0, nsim=1000, ci=0.95))$data[c("Median","Min","Max")])
sim_dist2_ppg_pl <- cbind(simGG(coxsimLinear(m4pl, b="dist2", Xj=1, Xl=0, nsim=1000, ci=0.95))$data[c("Median","Min","Max")])


# Romania1
sim_prominence_ppg_ro1 <- cbind(simGG(coxsimLinear(m3ro1, b="prominence", Xj=1, Xl=0, nsim=1000, ci=0.95))$data[c("Median","Min","Max")])
sim_MinorParty_ppg_ro1 <- cbind(simGG(coxsimLinear(m3ro1, b="MinorParty", Xj=1, Xl=0, nsim=1000, ci=0.95))$data[c("Median","Min","Max")])
sim_govmaj_ppg_ro1 <- cbind(simGG(coxsimLinear(m3ro1, b="govmaj", Xj=govmaj_upper, Xl=govmaj_lower, nsim=1000, ci=0.95))$data[c("Median","Min","Max")])
sim_collapse_ppg_ro1 <- cbind(simGG(coxsimLinear(m3ro1, b="collapse", Xj=1, Xl=0, nsim=1000, ci=0.95))$data[c("Median","Min","Max")])
sim_Expulsion_ppg_ro1 <- cbind(simGG(coxsimLinear(m3ro1, b="Expulsion", Xj=1, Xl=0, nsim=1000, ci=0.95))$data[c("Median","Min","Max")])
sim_relcoswitcher_ppg_ro1 <- cbind(simGG(coxsimLinear(m3ro1, b="relcoswitcher", Xj=relcoswitcher_upper, Xl=0, nsim=1000, ci=0.95))$data[c("Median","Min","Max")])
sim_preel_ppg_ro1 <- cbind(simGG(coxsimLinear(m3ro1, b="preel", Xj=1, Xl=0, nsim=1000, ci=0.95))$data[c("Median","Min","Max")])
sim_govformation_ppg_ro1 <- cbind(simGG(coxsimLinear(m3ro1, b="govformation", Xj=1, Xl=0, nsim=1000, ci=0.95))$data[c("Median","Min","Max")])
sim_dist2_ppg_ro1 <- cbind(simGG(coxsimLinear(m4ro1, b="dist2", Xj=1, Xl=0, nsim=1000, ci=0.95))$data[c("Median","Min","Max")])


# Romania2
sim_prominence_ppg_ro2 <- cbind(simGG(coxsimLinear(m3ro2, b="prominence", Xj=1, Xl=0, nsim=1000, ci=0.95))$data[c("Median","Min","Max")])
sim_MinorParty_ppg_ro2 <- cbind(simGG(coxsimLinear(m3ro2, b="MinorParty", Xj=1, Xl=0, nsim=1000, ci=0.95))$data[c("Median","Min","Max")])
sim_govmaj_ppg_ro2 <- cbind(simGG(coxsimLinear(m3ro2, b="govmaj", Xj=govmaj_upper, Xl=govmaj_lower, nsim=1000, ci=0.95))$data[c("Median","Min","Max")])
sim_collapse_ppg_ro2 <- cbind(simGG(coxsimLinear(m3ro2, b="collapse", Xj=1, Xl=0, nsim=1000, ci=0.95))$data[c("Median","Min","Max")])
sim_Expulsion_ppg_ro2 <- cbind(simGG(coxsimLinear(m3ro2, b="Expulsion", Xj=1, Xl=0, nsim=1000, ci=0.95))$data[c("Median","Min","Max")])
sim_relcoswitcher_ppg_ro2 <- cbind(simGG(coxsimLinear(m3ro2, b="relcoswitcher", Xj=relcoswitcher_upper, Xl=0, nsim=1000, ci=0.95))$data[c("Median","Min","Max")])
sim_preel_ppg_ro2 <- cbind(simGG(coxsimLinear(m3ro2, b="preel", Xj=1, Xl=0, nsim=1000, ci=0.95))$data[c("Median","Min","Max")])
sim_govformation_ppg_ro2 <- cbind(simGG(coxsimLinear(m3ro2, b="govformation", Xj=1, Xl=0, nsim=1000, ci=0.95))$data[c("Median","Min","Max")])
sim_dist2_ppg_ro2 <- cbind(simGG(coxsimLinear(m4ro2, b="dist2", Xj=1, Xl=0, nsim=1000, ci=0.95))$data[c("Median","Min","Max")])
##################################################################################



##################################################################################
# Create a data frame for all the simulation results
# Outcome variable: PPG entry
entry_plot_data <- rbind(
  # Prominence
  data.frame(Variable = "prominence", Country = "lt", sim_prominence_entry_lt),
  data.frame(Variable = "prominence", Country = "pl", sim_prominence_entry_pl),
  data.frame(Variable = "prominence", Country = "ro1", sim_prominence_entry_ro1),
  data.frame(Variable = "prominence", Country = "ro2", sim_prominence_entry_ro2),
  
  # MinorParty
  data.frame(Variable = "MinorParty", Country = "lt", sim_MinorParty_entry_lt),
  data.frame(Variable = "MinorParty", Country = "pl", sim_MinorParty_entry_pl),
  data.frame(Variable = "MinorParty", Country = "ro1", sim_MinorParty_entry_ro1),
  data.frame(Variable = "MinorParty", Country = "ro2", sim_MinorParty_entry_ro2),
  
  # dist2
  data.frame(Variable = "dist2", Country = "lt", sim_dist2_entry_lt),
  data.frame(Variable = "dist2", Country = "pl", sim_dist2_entry_pl),
  data.frame(Variable = "dist2", Country = "ro1", sim_dist2_entry_ro1),
  data.frame(Variable = "dist2", Country = "ro2", sim_dist2_entry_ro2),
  
  # govmaj
  data.frame(Variable = "govmaj", Country = "lt", sim_govmaj_entry_lt),
  data.frame(Variable = "govmaj", Country = "pl", sim_govmaj_entry_pl),
  data.frame(Variable = "govmaj", Country = "ro1", sim_govmaj_entry_ro1),
  data.frame(Variable = "govmaj", Country = "ro2", sim_govmaj_entry_ro2),
  
  # collapse
  data.frame(Variable = "collapse", Country = "lt", sim_collapse_entry_lt),
  data.frame(Variable = "collapse", Country = "pl", sim_collapse_entry_pl),
  data.frame(Variable = "collapse", Country = "ro1", sim_collapse_entry_ro1),
  data.frame(Variable = "collapse", Country = "ro2", sim_collapse_entry_ro2),
  
  # Expulsion
  data.frame(Variable = "Expulsion", Country = "lt", sim_Expulsion_entry_lt),
  data.frame(Variable = "Expulsion", Country = "pl", sim_Expulsion_entry_pl),
  data.frame(Variable = "Expulsion", Country = "ro1", sim_Expulsion_entry_ro1),
  data.frame(Variable = "Expulsion", Country = "ro2", sim_Expulsion_entry_ro2),
  
  # relcoswitcher
  data.frame(Variable = "relcoswitcher", Country = "lt", sim_relcoswitcher_entry_lt),
  data.frame(Variable = "relcoswitcher", Country = "pl", sim_relcoswitcher_entry_pl),
  data.frame(Variable = "relcoswitcher", Country = "ro1", sim_relcoswitcher_entry_ro1),
  data.frame(Variable = "relcoswitcher", Country = "ro2", sim_relcoswitcher_entry_ro2),
  
  # preel
  data.frame(Variable = "preel", Country = "lt", sim_preel_entry_lt),
  data.frame(Variable = "preel", Country = "pl", sim_preel_entry_pl),
  data.frame(Variable = "preel", Country = "ro1", sim_preel_entry_ro1),
  data.frame(Variable = "preel", Country = "ro2", sim_preel_entry_ro2),
  
  # govformation
  data.frame(Variable = "govformation", Country = "lt", sim_govformation_entry_lt),
  data.frame(Variable = "govformation", Country = "pl", sim_govformation_entry_pl),
  data.frame(Variable = "govformation", Country = "ro1", sim_govformation_entry_ro1),
  data.frame(Variable = "govformation", Country = "ro2", sim_govformation_entry_ro2)
)

# Add event type
entry_plot_data$EventType <- "Entry"

# Outcome variable: PPG creation
ppg_plot_data <- rbind(
  # Prominence
  data.frame(Variable = "prominence", Country = "lt", sim_prominence_ppg_lt),
  data.frame(Variable = "prominence", Country = "pl", sim_prominence_ppg_pl),
  data.frame(Variable = "prominence", Country = "ro1", sim_prominence_ppg_ro1),
  data.frame(Variable = "prominence", Country = "ro2", sim_prominence_ppg_ro2),
  
  # MinorParty
  data.frame(Variable = "MinorParty", Country = "lt", sim_MinorParty_ppg_lt),
  data.frame(Variable = "MinorParty", Country = "pl", sim_MinorParty_ppg_pl),
  data.frame(Variable = "MinorParty", Country = "ro1", sim_MinorParty_ppg_ro1),
  data.frame(Variable = "MinorParty", Country = "ro2", sim_MinorParty_ppg_ro2),
  
  # dist2
  data.frame(Variable = "dist2", Country = "lt", sim_dist2_ppg_lt),
  data.frame(Variable = "dist2", Country = "pl", sim_dist2_ppg_pl),
  data.frame(Variable = "dist2", Country = "ro1", sim_dist2_ppg_ro1),
  data.frame(Variable = "dist2", Country = "ro2", sim_dist2_ppg_ro2),
  
  # govmaj
  data.frame(Variable = "govmaj", Country = "lt", sim_govmaj_ppg_lt),
  data.frame(Variable = "govmaj", Country = "pl", sim_govmaj_ppg_pl),
  data.frame(Variable = "govmaj", Country = "ro1", sim_govmaj_ppg_ro1),
  data.frame(Variable = "govmaj", Country = "ro2", sim_govmaj_ppg_ro2),
  
  # collapse
  data.frame(Variable = "collapse", Country = "lt", sim_collapse_ppg_lt),
  data.frame(Variable = "collapse", Country = "pl", sim_collapse_ppg_pl),
  data.frame(Variable = "collapse", Country = "ro1", sim_collapse_ppg_ro1),
  data.frame(Variable = "collapse", Country = "ro2", sim_collapse_ppg_ro2),
  
  # Expulsion
  data.frame(Variable = "Expulsion", Country = "lt", sim_Expulsion_ppg_lt),
  data.frame(Variable = "Expulsion", Country = "pl", sim_Expulsion_ppg_pl),
  data.frame(Variable = "Expulsion", Country = "ro1", sim_Expulsion_ppg_ro1),
  data.frame(Variable = "Expulsion", Country = "ro2", sim_Expulsion_ppg_ro2),
  
  # relcoswitcher
  data.frame(Variable = "relcoswitcher", Country = "lt", sim_relcoswitcher_ppg_lt),
  data.frame(Variable = "relcoswitcher", Country = "pl", sim_relcoswitcher_ppg_pl),
  data.frame(Variable = "relcoswitcher", Country = "ro1", sim_relcoswitcher_ppg_ro1),
  data.frame(Variable = "relcoswitcher", Country = "ro2", sim_relcoswitcher_ppg_ro2),
  
  # preel
  data.frame(Variable = "preel", Country = "lt", sim_preel_ppg_lt),
  data.frame(Variable = "preel", Country = "pl", sim_preel_ppg_pl),
  data.frame(Variable = "preel", Country = "ro1", sim_preel_ppg_ro1),
  data.frame(Variable = "preel", Country = "ro2", sim_preel_ppg_ro2),
  
  # govformation
  data.frame(Variable = "govformation", Country = "lt", sim_govformation_ppg_lt),
  data.frame(Variable = "govformation", Country = "pl", sim_govformation_ppg_pl),
  data.frame(Variable = "govformation", Country = "ro1", sim_govformation_ppg_ro1),
  data.frame(Variable = "govformation", Country = "ro2", sim_govformation_ppg_ro2)
)

# Add event type
ppg_plot_data$EventType <- "PPG Creation"

# Combine all data
all_plot_data <- rbind(entry_plot_data, ppg_plot_data)

# Full word country labels
country_labels <- c(
  lt = "Lithuania",
  pl = "Poland",
  ro1 = "Romania1",
  ro2 = "Romania2"
)
all_plot_data$CountryLabel <- factor(country_labels[all_plot_data$Country], 
                                     levels = c("Romania2","Romania1","Poland","Lithuania"),
                                     labels = c("RomaniaMMP","RomaniaPR","Poland","Lithuania"))

# More informative variable labels
all_plot_data$VariableLabel <- factor(all_plot_data$Variable, 
                                      levels = c("prominence", "MinorParty","dist2",
                                                 "govmaj", "collapse", 
                                                 "Expulsion", "relcoswitcher", "preel","govformation"),
                                      labels = c("Personal Vote", "Minor Party", 
                                                 "Policy proximity",
                                                 "Gov. Majority", 
                                                 "Collapse", "Expulsion", "(Rel.) Co-Switcher", 
                                                 "Pre-Election","Gov. Formation"))


all_plot_data$Median <- log(all_plot_data$Median)
all_plot_data$Min <- log(all_plot_data$Min)
all_plot_data$Max <- log(all_plot_data$Max)
##################################################################################


##################################################################################
########################### Figure 4 ############################################
main_effects <- all_plot_data %>% filter(Variable %in% c("prominence", "MinorParty"))
main_effects$VariableLabel <- relevel(main_effects$VariableLabel,ref="Minor Party")
main_effects$hr <- as.character(round((exp(main_effects$Median)-1)*100,0))

p <- ggplot(main_effects, aes(x = VariableLabel, y = Median, color = CountryLabel)) +
  geom_point(position = position_dodge(width = 0.5), size = 3) +
  geom_errorbar(aes(ymin = Min, ymax = Max), 
                position = position_dodge(width = 0.5), width = 0.3) +
  geom_text(aes(label = hr, y = Median), 
            position = position_dodge(width = 0.5), 
            size = 4,
            vjust = -0.8) +
  facet_wrap(~ EventType, ncol = 1) +
  coord_flip() +
  labs(
    title = "",
    x = "Variable",
    y = "Hazard Ratio (log scale)",
    color = "Country",
    size=5
  ) +
  scale_color_grey(
    breaks = c("Lithuania", "Poland", "RomaniaPR", "RomaniaMMP")
  ) +
  theme_bw() +
  theme(
    legend.position = "bottom",
    legend.text = element_text(size=13),
    panel.grid.major = element_line(color = "gray90"),
    panel.grid.minor = element_line(color = "gray95"),
    strip.text = element_text(face = "bold", size = 12),
    axis.text.y = element_text(size = 12)
  ) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "darkgray", alpha = 0.7) #+
ggsave("figure_4.png", height = 6, width = 7, dpi = 700)
##################################################################################


##################################################################################
########################### Figure 6 ############################################
controls <- all_plot_data %>% filter(!(Variable %in% c("prominence","MinorParty")))
controls$CountryLabel <- factor(country_labels[controls$Country], 
                                     levels = c("Lithuania","Poland","Romania1","Romania2"),
                                     labels = c("LT","PL","RO-PR","RO-MMP"))
controls$hr <- as.character(round((exp(controls$Median)-1)*100,0))

# Add a new column to determine statistical significance
controls <- controls %>%
  mutate(Significant = !(Min <= 0 & Max >= 0))

ggplot(controls, aes(x = CountryLabel, y = Median, color = Significant)) +
  geom_point(size = 3) +
  geom_errorbar(aes(ymin = Min, ymax = Max), width = 0.25) +
  geom_text(
    data = subset(controls, Significant == TRUE),
    aes(label = format(hr, digits = 2, nsmall = 2)),
    vjust = -1,
    size = 3
  ) +
  facet_grid(EventType ~ VariableLabel, scales = "free_y") +
  labs(
    title = "",
    x = "",
    y = "Hazard Ratio (log scale)"
  ) +
  scale_color_manual(values = c("gray", "black")) +
  theme_bw() +
  theme(
    legend.position = "none",
    axis.text.x = element_text(angle = 45, hjust = 1),
    panel.grid.major = element_line(color = "gray90"),
    panel.grid.minor = element_line(color = "gray95"),
    strip.text = element_text(size = 9),
    axis.title.x = element_text(size = 10)
  ) +
  scale_y_continuous(labels = scales::number_format(accuracy = 0.01)) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "darkgray", alpha = 0.7)
ggsave("figure_6.png", height = 8, width = 12, dpi = 900)
##################################################################################



##################################################################################
###################### FIGURE 5 ##################################################

# Joint Cox model
m1j <- coxph(list(Surv(start,stop,event,origin=0) ~ prominence + MinorParty + govmaj + collapse + Expulsion + preel+ govformation, 1:2 ~ coswitcher, 1:3 ~ relcoswitcher),pl,id=spellid,cluster=spellid, robust=T)

# function to get cumulative incidence functions
f_inc <- function(model,spells,spellid,start,stop){
  tmp <- data.frame(start=start,stop=stop,event=0,new[new$spellid==spellid & new$start==start,])
  res <-survfit(m1j,newdata= tmp,se.fit=T) 
  figtime<-rep(res$time,2)
  figps1<-res$pstate[,,2]
  figps2<-res$pstate[,,3]
  fallp<-c(figps1,figps2)
  fdata<-data.frame(figtime,fallp,Event=c(rep("PPG entry",length(figtime)/2),rep("PPG creation",length(figtime)/2)))
  titleinfo <- paste(tmp$country[1],tmp$name[1],
                    "Personal Vote =",tmp$prominence[1],
                     "Minor Party =",tmp$MinorParty[1])
  tmp$collapse[tmp$collapse==0] <- "NO"
  tmp$collapse[tmp$collapse==1] <- "YES"
  tmp$MinorParty[tmp$MinorParty==0] <- "NO"
  tmp$MinorParty[tmp$MinorParty==1] <- "YES"
  plot<-ggplot(data=fdata,aes(x=figtime,y=fallp,color=Event,linetype=Event))+
    geom_step(direction="hv",size=1)+theme_bw()+labs(x="Days", y="",title=titleinfo)+
    xlim(0,1400)+ylim(0,1)+
    theme(legend.position="none",plot.title = element_text(size=8))+
    scale_color_manual(values=c("gray","black"))+
    scale_linetype_manual(values=c(2,1))
  if (tmp$name[1]%in%c("Rucinski Andrzej")){
    plot <- plot + geom_vline(aes(xintercept=tmp$lastday[1]),linetype=1,col="black")
  }
  if (tmp$name[1]%in%c("Zalek Jacek","Janowski Gabriel")){
    plot <- plot + geom_vline(aes(xintercept=tmp$lastday[1]),linetype=2,col="gray")
  }
  data <- tmp[c("name","lastday", "coswitcher","relcoswitcher", "collapse","MinorParty","govmaj")]
  data$event <- as.character(spells$event[spells$spellid==tmp$spellid[1]])
  data$event[data$event=="0"] <- "None"
  data$event[data$event=="1"] <- "PPG entry"
  data$event[data$event=="2"] <- "PPG creation"
  return(list(plot,data))
}
plot1pl <- f_inc(m1j,spells,"PL Rucinski Andrzej 2007-04-26",29,30)[[1]]
plot2pl <- f_inc(m1j,spells,"PL Zalek Jacek 2014-07-14",29,30)[[1]]
plot3pl <- f_inc(m1j,spells,"PL Dorn Ludwik 2011-10-09",29,30)[[1]]
plot4pl <- f_inc(m1j,spells,"PL Janowski Gabriel 2005-06-30",29,30)[[1]]

incidence <- ggarrange(plot1pl,plot3pl,plot2pl,plot4pl,
                       common.legend = TRUE, legend="bottom",ncol=2,nrow=2) +
  theme_bw(base_family = "serif")
ggsave("figure_5.png", height = 6, width = 8, dpi = 800)
##############################################################################################






##############################################################################################
##############################################################################################
##############################################################################################
#########################APPENDICES##########################################################
## Appendix 2: Descriptive statistics for explanatory and control variables
ltdata1 <- new %>% filter(duplicated(spellid)==FALSE) %>%
  rename("PersonalVote"="prominence","MinorParty"="MinorParty",
         "PPGcollapse"="collapse","Expulsion"="Expulsion",
         "CoSwitcher"="coswitcher","RelCoSwitcher"="relcoswitcher") %>%
  filter(country=="LT")
pldata1 <- new %>% filter(duplicated(spellid)==FALSE) %>%
  rename("PersonalVote"="prominence","MinorParty"="MinorParty",
         "PPGcollapse"="collapse","Expulsion"="Expulsion",
         "CoSwitcher"="coswitcher","RelCoSwitcher"="relcoswitcher") %>%
  filter(country=="PL")
rodata11 <- new %>% filter(duplicated(spellid)==FALSE) %>%
  rename("PersonalVote"="prominence","MinorParty"="MinorParty",
         "PPGcollapse"="collapse","Expulsion"="Expulsion",
         "CoSwitcher"="coswitcher","RelCoSwitcher"="relcoswitcher") %>%
  filter(country=="RO1")
rodata12 <- new %>% filter(duplicated(spellid)==FALSE) %>%
  rename("PersonalVote"="prominence","MinorParty"="MinorParty",
         "PPGcollapse"="collapse","Expulsion"="Expulsion",
         "CoSwitcher"="coswitcher","RelCoSwitcher"="relcoswitcher") %>%
  filter(country=="RO2")
descdata1 <- rbind(ltdata1,pldata1,rodata11,rodata12) 

datasummary((PersonalVote+MinorParty+PPGcollapse+
               Expulsion+CoSwitcher+RelCoSwitcher) ~ 
              (Min + Mean + Median + SD + Max + N), data=ltdata1,
            title= "Descriptive statistics: time-invariant variables: Lithuania",
            output="kableExtra") %>% kable_styling(full_width = F, latex_options = "hold_position")
datasummary((PersonalVote+MinorParty+PPGcollapse+
               Expulsion+CoSwitcher+RelCoSwitcher) ~ 
              (Min + Mean + Median + SD + Max + N), data=pldata1,
            title= "Descriptive statistics: time-invariant variables: Poland",
            output="kableExtra") %>% kable_styling(full_width = F, latex_options = "hold_position")
datasummary((PersonalVote+MinorParty+PPGcollapse+
               Expulsion+CoSwitcher+RelCoSwitcher) ~ 
              (Min + Mean + Median + SD + Max + N), data=rodata11,
            title= "Descriptive statistics: time-invariant variables: Romania-PR",
            output="kableExtra") %>% kable_styling(full_width = F, latex_options = "hold_position")
datasummary((PersonalVote+MinorParty+PPGcollapse+
               Expulsion+CoSwitcher+RelCoSwitcher) ~ 
              (Min + Mean + Median + SD + Max + N), data=rodata12,
            title= "Descriptive statistics: time-invariant variables: Romania-MMP",
            output="kableExtra") %>% kable_styling(full_width = F, latex_options = "hold_position")



## Appendix 4: Robustness checks

# alternative measures of personal vote
# without distance variable, ENTRY
m1lt_1 <- coxph(Surv(start,stop,entryevent) ~ logprefshare + smdsurplus + MinorParty + govmaj + collapse + Expulsion + coswitcher + preel+ govformation,lt,id=spellid,cluster=spellid, robust=T)
m1pl_1 <- coxph(Surv(start,stop,entryevent) ~ logprefshare + MinorParty + govmaj + collapse + Expulsion + coswitcher + preel+ govformation,pl,id=spellid,cluster=spellid, robust=T)
m1ro2_1 <- coxph(Surv(start,stop,entryevent) ~ smdsurplus + MinorParty + govmaj + collapse + Expulsion + coswitcher + preel+ govformation,ro2,id=spellid,cluster=spellid, robust=T)
m3lt_1 <- coxph(Surv(start,stop,PPGcrevent) ~ logprefshare + smdsurplus + MinorParty + govmaj + collapse + Expulsion + coswitcher + preel+ govformation,lt,id=spellid,cluster=spellid, robust=T)
m3pl_1 <- coxph(Surv(start,stop,PPGcrevent) ~ logprefshare + MinorParty + govmaj + collapse + Expulsion + coswitcher + preel+ govformation,pl,id=spellid,cluster=spellid, robust=T)
m3ro2_1 <- coxph(Surv(start,stop,PPGcrevent) ~ smdsurplus + MinorParty + govmaj + collapse + Expulsion + coswitcher + preel+ govformation,ro2,id=spellid,cluster=spellid, robust=T)


# fixed effects for parliamentary terms
m2lt <- update(m1lt,.~.+dist2+country_term)
m2pl <- update(m2lt,data=pl)
m2pl <- update(m2pl,.~.+tMinorParty+tgovmaj+tpreel+country_term)
m2ro1 <- update(m2lt,data=ro1)
m2ro1 <- update(m2ro1,.~.+tgovmaj+tcollapse+trelcoswitcher+country_term)
m2ro2 <- update(m2lt,data=ro2)
m4lt <- update(m3lt,.~.+dist2+country_term)
m4pl <- update(m4lt,data=pl)
m4ro1 <- update(m4lt,data=ro1)
m4ro2 <- update(m4lt,data=ro2)
##############################################################################################
##############################################################################################
##############################################################################################




